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Abstract. We derive a set of algorithms for simulating the diffusion- 
limited growth of faceted crystals using local cellular automata. This 
technique has been shown to work well in reproducing realistic crystal 
£o ! morphologies, and the present work provides a more rigorous phys- 

ical foundation that connects the numerical code to the physics of 
attachment kinetics and diffusion dynamics. We then apply these al- 
gorithms to examine a novel morphological transition in the growth 
of thin plate-like crystals. 

Y ; 1 Introduction 

£h ' The formation of complex structures during solidification often results from 

a subtle interplay of nonequilibrium, nonlinear processes, for which seemingly 
small changes in molecular dynamics at the nanoscale can produce profound 
morphological changes at all scales. One popular example of this phenomenon 
is the formation of snow crystals, which are ice crystals that grow from wa- 
ter vapor in an inert background gas. Although this is a relatively simple, 
monomolecular system, snow crystals display a remarkable variety of columnar 

\^ . and plate-like forms, and much of the phenomenology of their growth remains 

C^l ' poorly understood, even at a qualitative level [1] . 

r — . Viewed broadly, snow crystal structures result from diffusion-limited crystal 

f^ ' growth in the presence of strongly anisotropic attachment kinetics, and similar 

circumstances occur in a large class of solidification problems. We can break 

^-^ ' such problems down into two main physical components: 1) the attachment 

kinetics that describes the molecular growth dynamics at the solid surface, and 
2) particle transport via diffusion to the growing crystal. In the present paper 
we focus on the latter problem - numerically solving the diffusion equation to 
model faceted crystal growth. 

Computational models of diffusion-limited solidification have typically been 
divided into two broad camps - 'front-tracking' models, in which one keeps track 
of the solidification interface explicitly (e.g. [21 02 0]), and 'phase- field' mod- 
els, in which the solidification front is numerically smoothed and not explicitly 
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tracked (e.g. [SJ[5]). To date, both techniques have had considerable success 
in modeling simple dendritic growth, but neither has been able to satisfacto- 
rily model complex morphological structures in the presence of strong faceting 
[51 M QH] > owing to the appearance of dynamical and numerical instabilities. 

A third computation technique - local cellular automata (LCA) - has recently 
been applied to the problem of modeling faceted crystal growth with excellent 
initial success, and LCA models have produced realistic-looking snow crystal 
growth simulations in both 2D [TTJ[T2] and 3D [T3]. To date, however these local 
lattice models have incorporated largely ad hoc growth rules (albeit physically 
motivated to some extent), so the connection between the resulting simulations 
and real crystal growth remains tenuous. These algorithms yielded structures 
that bear a striking resemblance to many features seen in real snow crystals, 
but the similarities are largely qualitative. As a result, these models cannot be 
directly compared with crystal growth measurements to provide quantitative 
insights into growth dynamics. 

The goal of this paper is to produce more accurate, physically derived rules 
for modeling with cellular automata, thus providing a direct connection between 
the numerical code and the physics of attachment kinetics and diffusion dynam- 
ics governing faceted crystal growth. The algorithms derived below provide a 
more rigorous physical foundation for using LCA techniques in crystal growth 
models. After deriving growth algorithms below for the case of snow crystal 
formation, we then examine the properties of the LCA model and describe a 
novel morphological transition seen in plate-like growth. 

2 Crystal Growth Dynamics 

We focus our discussion on the growth of snow crystals, as this is perhaps 
the most studied example of strongly faceted, diffusion-limited crystal growth. 
Having such a focus allows a more direct look at the relevant physics, which 
becomes somewhat obscured in a more general discussion. Extending the al- 
gorithms to other materials should be relatively straightforward, provided the 
growth dynamics are of a similar character. 

We begin by assuming the " standard model" of snow crystal growth [1] , for 
which we can first write the growth velocity normal to the surface in terms of 
the Hcrtz-Knudscn formula 



Csat kT 

v n = a \ a sur f (1) 
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where the latter expression defines the velocity Vkin- In this expression kT is 
Boltzmann's constant times temperature, m is the mass of a water molecule, 
CsoUd = Pice/m is the number density for ice, a surf = (c surf - c sa t)/c sa t is the 
supersaturation just above the growing surface, c sur f is the water vapor number 
density at the surface, and c sat (T) is the equilibrium number density above a flat 



ice surface. The dimensionless parameter a is known as the condensation coeffi- 
cient, and it embodies the surface physics that governs how water molecules are 
incorporated into the ice lattice, collectively known as the attachment kinetics. 

The attachment kinetics can be nontrivial, so in general a will depend on 
T, a sur f, and perhaps on the surface structure and geometry, surface chemistry, 
and other factors. If molecules striking the surface are instantly incorporated 
into it, then a = 1; otherwise we must have a < 1. The appearance of crystal 
facets indicates that the growth is limited in part by attachment kinetics, so we 
must have a < 1 on faceted surfaces. 

This model assumes that the attachment kinetics are purely local, without 
significant large-scale surface diffusion, especially around corners between facets. 
It also assumes that all complex aspects of the molecular dynamics at the ice 
surface, which may include surface melting and any number of other details, can 
be absorbed into some parameterization of a as a function of external parame- 
ters. There is currently no evidence that this "standard model" is incorrect pQ, 
but at the same time we cannot prove that the attachment kinetics is always 
well described by Equation [TJ For the present discussion, we will assume that 
this expression is valid for circumstances relevant to snow crystal growth. 

Particle transport through the air surrounding a growing crystal is described 
by the diffusion equation 

g = DV*c (2, 

where c(x) is the water molecule number density surrounding the crystal and D 
is the diffusion constant. The timescale for diffusion to adjust the vapor concen- 
tration in the vicinity of a crystal is Tdif fusion ~ R 2 /D, where R is a characteris- 
tic crystal size. This is to be compared with the growth time, T gro wth ~ 2R/v n , 
where v n is the growth velocity of the solidification front normal to the sur- 
face. The ratio of these two timescales is the Peclet number, p = Rv n /2D. 
For typical growth rates of snow crystals we find p < 10 -5 , which means that 
diffusion adjusts the particle density around the crystal much faster than the 
crystal shape changes. In this case the diffusion equation reduces to Laplace's 
equation, V 2 c = 0, which must be solved with the appropriate boundary con- 
ditions. Using this slow-growth limit often simplifies the problem considerably 
in comparison to much of the literature on diffusion- limited solidification [14] . 
The continuity equation at the interface gives 
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where a(x) = [c(x)—c sa t}/c S at and we are assuming the isothermal case, so c sa t is 
independent of spatial position. The latter assumption means we will be ignoring 
effects that arise when a growing crystal experiences a temperature increase 
associated with the latent heat of solidification. These effects are generally small 
and produce results that are similar to a simple decrease in the supersaturation 
surrounding the crystal pQ. Thus we believe we will not be sacrificing much 
interesting physics by using the isothermal approximation. We then write the 



diffusion equation in terms of the supersaturation field (since c sa t is constant) 
as 

r)rr 

DV 2 a (4) 
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The attachment coefficient a is not well known for ice, but on facet surfaces 
it appears to be well described by nucleation-limited growth [TJ [T51 [TB] , which 
gives 

a(cr) ~ Aexp (— <Jo/cr) (5) 

where A and do are parameters that depend on temperature and are different 
for the basal and prism facets. For growth governed by spiral dislocations, we 
expect [T7] 

a(a) « Co (6) 

where C is independent of a (but depends on external parameters such as tem- 
perature), which gives v n ~ cr 2 . 

3 Derivation of Growth Algorithms 

3.1 The ID Problem 

Because of its simplicity, it is instructive to begin with the ID crystal growth 
problem, for which we assume a linear array of grid elements, or pixels, of width 
Ax. Each pixel is assumed to be filled with either ice or air at any given time, 
and each air pixel has a well-defined supersaturation cr(x,t) associated with it. 
The discrete form of the ID diffusion equation is 



u(x, t + fit) — er(x, t) 
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a(x + dx) — 2a(x) + o(x — dx) 



dx 2 



dt 

which yields 

a(x, t + At) = Ato(x - Ax, t) + (1 - 2Ar) a(x, t) + Ato(x + Ax, t) (8) 

where At is the physical time step in our simulation code and Ax is the grid 
spacing (assumed uniform). The quantity 

DAt 

Ar = 7^ (9) 

(Ax) 

is a dimensionless time-step parameter. We note that (Ax) / D is one natu- 
ral time scale of the problem, equal the time required for diffusion to adjust 
the supersaturation field over a distance Ax. We define r = J2 ^ T to be the 
dimensionless physical time in our problem. 

We use Equation [8] to effectively relax the supersaturation field to a suit- 
able solution of the diffusion equation. Choosing Ar too small will require an 



excessive amount of computer time, while taking At > 1 leads to instabilities 
near sharp edges (e.g., at physical boundaries) in a. It appears that At = 0.5 
is close to optimal, yielding the simple algorithm 

a(x, t + At) = ~a(x - Ax, r) + -a(x + Ax, r) (10) 

for the ID case. This expression is used to propagate a in time for regions 
away from the growing crystal surface, where diffusion alone affects the particle 
dynamics. 

Our model crystal grows as water molecules diffuse in from the outer bound- 
aries of the system and attach to the ice surface. For boundary pixels (i.e., air 
pixels that are adjacent to ice pixels), we see physically that the ice surface ef- 
fectively "drains" the supersaturation as water vapor condenses, which in turn 
causes the solidification front to move with a velocity v = avki n a. This "drain- 
ing" reduces the supersaturation of boundary pixels according to 

g(t + At) = cr(l - aA£,Ar) (11) 

where Af = Ax/ 'X is a dimensionless pixel size in the problem and 

Csot D 
X = — (12) 
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is the natural physical scale of the problem. Here D a i r = 2x 10~ 5 m/sec 2 is the 
diffusion constant for air at one atmospheric pressure, and the numerical value 
pertains to ice crystal growth. 

We can combine Equation II II with Equation [5] to create a new propagation 
algorithm for boundary pixels 

a(x, t + At) = Ara so ud (r) + (1 - 2 At) a(x, t) + At<t(x + Ax,t) (13) 

where the boundary pixel is at x and the ice pixel is at x = x — Ax, and 

<r a oiid(T) = cr(x, r)(l - aA£_) (14) 

With this, the same algorithm can be used for all air pixels, including boundary 
pixels, provided one substitutes <J S olid m the appropriate place for boundary 
pixels. 

For the conversion of boundary pixels to ice, mass conservation suggests that 
we define an accumulated mass parameter A that begins as A = for each air 
pixel. For boundary pixels, we increment A using 

A^A + acrAA (15) 



for each time step, where 

AA = _2f^lA T A£ (16) 

C-solid 

When A > 1, that pixel converts from air to ice. 

For ice we have c sa t/c so ud w 10~ 6 , while the dimensionless parameters a, 
At, a, A£ will all be of order unity or perhaps substantially smaller. Thus it 
would take more than a million time steps before a boundary pixel becomes 
an ice pixel. This situation reflects the small Peclet number in our problem, 
so the growth is much slower than the time it takes diffusion to adjust the 
supersaturation field. We speed up the code by taking 

AA = AAtA£ (f 7) 

where A is an input constant, discussed further below. 

3.2 The 2D problem in Cartesian Coordinates 

Following the ID analysis, the discrete form of the 2D diffusion equation gives 

<t(t+At) = At [a(x - Ax) + <j(x + Ax) + <j(y - Ay) + cr(y + Ay)] + (1 - 4Ar) a 

(18) 
where we have dropped indices when doing so does not confuse the expression. 
If we take our time parameter to be At = 1/4, this becomes 

(t(t + At) = - [a(x - Ax) + o(x + Ax) + a(y - Ay) + <r(y + Ay)} (19) 
For boundary pixels, we again use this expression and substitute 

(Jsolid = <t(1 - aAO (20) 

for each neighboring ice pixel, as with the ID case, and we have further assumed 
a uniform grid with Ax = Ay. 

We again define an accumulated mass parameter A and increment it with 

A^A + ^cmtAA (21) 

for each time step, where the sum is over the number of ice neighbors. When 
A > 1, that pixel converts from air to ice. 

In both these expressions we must choose a with some care, as its value will 
depend critically on the number and orientation of neighboring solid pixels. We 
label a boundary pixel with (N x ,N y ), where N x is the number of neighboring 
ice pixels in the x direction and N y is the number of ice neighbors in the y 
direction. Both N x and N y can take values 0, 1, or 2, giving nine cases for 
(N x ,N y ). The cases are: 

(0,0) - the pixel is an air pixel 

(0,1) - one ice neighbor in the y direction, so a — a y , the physical value 
appropriate for a y facet surface. 



(1.0) - one ice neighbor in the x direction, so a — a x for an x facet surface. 

(1.1) - a kink site, where the growth will not be nucleation-limited, since 
the corner provides a source of molecular steps. We do not know a priori what 
value to use for a on this site, but assume a constant a = a\\. 

(0,2), (1,2), (2,0), (2,1), (2,2) - these are all unusual cases where the growth 
will be fast, so it shouldn't matter much what we choose for a, as long as it is 
large. 

We can index these possibilities with a single number by computing a bound- 
ary parameter B — 2N^ + Ny, where N x is the number of x neighbors (0, 1, or 
2) and N y is the number of y neighbors. We then have B — for an air pixel, 
B = 1 for a y facet, B = 2 for an x facet, B = 3 for a (1,1) kink location, and 
B > 3 for all other cases. 

If we consider the special case where a is equal to some constant value, 
independent of the orientation of the surface with respect to the crystal lattice, 
then the growth velocity should equal v — avkma for all surfaces. For the (01) 
or (10) facet surfaces 13 in this constant-a case, we take a x = a y = a, while an 
analysis of the growth of a (11) surface shows that we must take an = a/^/2 if 
the above algorithm is to produce the correct growth velocity. 

3.3 The 2D Problem in Cylindrical Coordinates 

The discrete version of the Laplacian in cylindrical coordinates yields 



ct(t+At) = At 



Ar\ / Ar 

1 a(r- Ar) + ( H )a(r + Ar) + aiz - Az) + aiz + Az) +(1 - 4At) a 

2r J \ 2r 



(22) 
where we have assumed da/ 89 = to reduce the problem to 2D, thus yielding 
cylindrically symmetric crystal structures. The problem then proceeds essen- 
tially as with the 2D case in rectangular coordinates, and similar (1 — Ar/2r) 
correction factors are needed for o~ S oUd and AA. With the preferred value of 
At = 1/4 we have 



ct(t+At = - 
4 



Ar \ / Ar 

1 air - Ar) + H ) a(r + Ar) + aiz - Az) + aiz + Az) 

2r I \ 2r 



. (23) 
In our model, we define the first row of pixels to have r = 0, for which we 

cannot evaluate the above expression. Instead we revert to the Cartesian form 

of the Laplacian to give 

air + At) = At [4cr(r + Ar) + a(z - Az) + a(z + Az)} + (1 - 6At) a (24) 

for those pixels. If we take Ar = 1/4 (the same as for r ^ 0), this becomes 

o{t + At) = - [Aa(r + Ar) + a(z - Az) + a(z + Az)} a (25) 



2 The notation here - integers in parentheses without commas - refers to the usual Miller 
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notation - integers in parentheses with commas - which we used to label boundary pixels. 



and the negative term leads to instabilities. We could correct this by choosing 
a smaller Ar, but this would make the code run more slowly. An alternative is 
to use Ar = 1/6 for the r = special case only, giving 

<j(t + At) = - [4cr(r + Ar) + a(z ~ Az) + a(z + Az)} (26) 

for r — pixels. This isn't a perfect solution, but we find it does not cause 
significant problems in the code, because it is only applied to a single row of 
pixels. 

3.4 The 2D Problem in Hexagonal Coordinates 

We can model planar snow crystal growth using a hexagonal lattice, which can 
be mapped onto a rectangular lattice as shown in Figure 1 [13]. The discrete 
diffusion equation becomes 

2 6 

<t(t + At) = -Ar^V t + (l-4Ar)cr (27) 

O 

where the index i refers to the six nearest neighbors around each pixel. For 
boundary pixels, we again use this expression and substitute 

VsolU = <r(l - «A£) (28) 

for each neighboring boundary pixel, as above, and we have taken Ax to be the 
distance between nearest neighbors. 

We again define an accumulated mass parameter A and increment it with 

A ""*" A + q Yl aaAX ( 29 ) 

o 

for each time step, where the sum is over the number of ice neighbors. When 
A > 1, that pixel converts from air to ice. 

We label each boundary pixel with (N), where N is the number of nearest 
neighbors: 

(0) - an air pixel 

(1) - a boundary pixel at the tip of a hexagon. The growth will be low here, 
and it likely doesn't matter much what value is chosen for a. The growth should 
be sufficient, however, to allow the crystal to grow from its initial seed. 

(2) - this case refers to the normal growth of a prism facet, and an anal- 
ysis shows that we much choose a = (v3/2) a pr ism so the growth equals 
v — a P rismVkin(J for this surface. 

(3) - this refers to growth at a kink site, and the choice of a here has a strong 
effect on the transition from faceted to branched growth [T3] . 

(4), (5), (6) - these are unusual cases where the growth is not much affected 
by the choice of a as long as it is large (of order unity) . 



3.5 The 3D Problem in Hexagonal Coordinates 

The discrete diffusion equation becomes 

2 6 2 

ct(t + At) = - At TV* + At TV,,- + (1 - 6Ar) a (30) 

O 

where the index i refers to the six nearest neighbors around each pixel in the 
basal (xy) plane and j refers to the two neighbors perpendicular to this plane 
(the z direction). For boundary pixels, we again use this expression and substi- 
tute 

v solid = c(l - aA() (31) 

for each neighboring boundary pixel, as above. We have taken Ax to be the 
distance between nearest neighbors, which is the same in the basal plane as in 
the z direction. 

We again define an accumulated mass parameter A and increment it with 

2 6 2 

A ~* A + q 5Z aaAX + Y^ aaAX (32) 

i=i j=i 

for each time step, where the sums are over ice neighbors. When A > 1, that 
pixel converts from air to ice. 

We label each boundary pixel with (Ni,Nj), where Ni is the number of 
nearest neighbors in the basal plane and Nj is the number of neighbors in the 
z direction. We will assume that neighbors in the basal plane are contiguous. 
Non-contiguous cases are unusual, and it shouldn't matter much how we assign 
growth rates in those cases, as long as a is large. 

(0,0) - an air pixel 

(1,0) - a boundary pixel at the tip of a hexagon, with no z neighbors. This 
is similar to the case (1) in the 2D hexagonal problem above. 

(2,0) - similar to the (2) case in the 2D hexagonal case above, and we take 

a = (vo/2) Uprism- 

(3.0) - growth at a kink site, similar to the (3) case above, and again the 
value of a chosen will have a strong effect on the transition from faceted to 
branched growth [13]. 

(0,1) - growth of the basal facet; so we take a = otbasai- 

(1.1) - not well determined, but faster than (1,0) or (0,1) 
(2,1) - not well determined, but faster than (1,1) or (2,0) 

(3,1) - again not well determined, but faster than (2,1). The values of a 
used for the (JVj, 1) sites will strongly affect the transition from faceted growth 
to structure formation on the basal facet [T3] . 

The values of a used for the remaining sites should be high, and the details 
will likely not affect the growth dynamics significantly. 



4 Using the Algorithms 

4.1 Intrinsic Anisotropics 

To test our code, we modeled the growth of spherical crystals with isotropic 
attachment kinetics, where there is a simple analytic result for the growth rate. 
We limited the growth in our models to small changes in radius, to minimize 
non-spherical growth that eventually arises from the Mullins-Sekerka instability 
[17] . Our numerical model followed the prescriptions above for the 2D problem 
in cylindrical coordinates, with a x = a y = a = constant for the facet surfaces 
and an = a/y/2 for a (1,1) kink site. All higher-order sites with B > 3 are 
irrelevant for this problem, because of the convex geometry of the spherical 
surface. 

We found that our code generated growth rates that were always a few 
percent larger than the analytic theory for all (reasonable) choices of a, A£, 
A, and other input parameters. Upon closer investigation, we found that our 
algorithm produces growth of a (12) surface that is approximately 8 percent 
faster than the (01), (10), or (11) surfaces (the latter surfaces all grow at the 
same rate, which, by design, is the correct rate). 

This demonstrates that our choice of a fixed rectangular grid, along with only 
a small number of growth rules, yields an intrinsic anisotropy in the growth al- 
gorithm. This anisotropy could be reduced by adding higher-order corrections 
that specify slightly different a values for different boundary pixels, depend- 
ing on the surface configuration beyond just the nearest neighbors, but such 
corrections would be difficult to implement. 

This intrinsic anisotropy should not be a significant effect for strongly faceted 
growth, but we expect that our LCA approach would not produce accurate 
morphological results for cases where the anisotropy in the attachment kinetics 
is less than approximately 10 percent. 

4.2 Scaling Behavior 

If we run a growth code and produce some complex crystal shape, the inter- 
pretation of our result still contains an ambiguity. The crystal size is given in 
pixels, where Aa; = A£Ao is the pixel size. The parameter A£ was fixed in 
the code, but X$ depends on the diffusion constant D, which is not otherwise 
specified. Similarly, a single time step in the code corresponds to a physical 
time 

At = ^r-Ar (33) 
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A 2 A^ 2 Ai 
D 
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Thus we see that the growth behavior at different air pressures (different D) 
is determined once we know the growth at a single pressure (provided a^ is the 
same at the different pressures). If the air pressure is half an atmosphere, the 
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growth morphology (however complex) will be the same as at one atmosphere, 
except in the former case the crystal will be double the size in double the time. 
This scaling behavior nicely explains why snow crystal morphology is generally 
simpler for smaller crystals and/or for lower air pressures, which has long been 
observed pQ. To my knowledge, this scaling behavior has not been identified in 
previous investigations of snow crystal growth. 

4.3 Limitations on the Grid Size 

The physical size of the grid is Ax — A£Ao, and there are limits to how coarse 
the grid can be without affecting the growth or causing instabilities in the code. 
Taking A£ > 1/a would cause <J so ud to become negative, which causes some 
concern in that it may produce instabilities in the code. With this limitation, 
the grid spacing could not be larger than Ax = X$/a. For air at a pressure of 
one atmosphere and a ~ 1, this gives Ax = 0.15 /im. Modeling a 1.5 mm snow 
crystal would then require a grid with at least 10,000 pixels on a side, which is 
something of a computational challenge. 

We can do better by realizing that a negative a so iid is not itself sufficient to 
cause instabilities. A closer look at the algorithm reveals that problems only 
begin when o{t + At) becomes negative in a single timestep, which happens 
when the supersaturation is " drained" to a negative value according to Equation 
ITT1 This puts a limitation of A£ > 1/qAt on the grid size, so we are able to 
use a coarser grid spacing if we also use a finer time step. Whatever the limit, 
it is important to note that our choice of grid spacing is not simply limited by 
a desire to produce a physically accurate model of crystal growth, but also by 
intrinsic instabilities in the code. 

Physically, we can gain some insights into these limitations from dendrite 
growth theory pQ. We have Xq « Rkin (the latter from Equation 28 in pQ), and 
a growing dendrite has a tip radius 

Rtip = Rkin (34) 

as 
Xo 
as 

where s is the dimensionless solvability parameter, which is of order unity for 
ice crystal growth pQ. The stability of the code (with At rs 1) thus limits 
the grid spacing to be no greater than the tip radius of a growing dendritic 
structure. Interestingly, it appears that the code can only function properly 
when the grid spacing is fine enough to allow the growth of physically realistic 
dendritic structures, the scale of which is given by solvability theory. 

4.4 Adaptive Time Steps 

We wish to choose A in Equation[17]as large as possible, so the code runs quickly, 
while not altering the growth appreciably. Our first criterion is that the growth 
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be slow over a single time step, so that it takes at least Nq ~ 10 time steps 
before a boundary pixel turns to ice, which means we must take 

A < ^A^JVo (35) 

Another criterion is that the Peclet number should be small, as dictated by the 
physics of the growth problem, for which we find 

/Ax\ At 1 

a< U)a?- (36) 

where R is a characteristic size of the crystal. Since Ri = R/Ax is typically 
larger than Nq, the latter requirement is the more stringent of the two. 
We can speed up the code further by taking 

A = A l \ (37) 

where A < 1 is a constant, -Ri. max is the current maximum size of the crystal 
(in pixels), and (ctc) max is the maximum product of a and a over all current 
boundary pixels. This speeds up the code considerably when the growth is 
strongly diffusion limited (so a <C Coo) or when a is small on the crystal surface. 
This choice of A means essentially using an adaptive time step, where the 
physical time for each step is equal to 

At = AA£ 2 ArAi (38) 

where 

27TTO —Cnlid 

Ato = D^± 39 

kl Csat 



D 



x (1 msec) 



Da 

and the latter value is for growth at T = — 15 C pQ. 

5 A Morphological Transition in Plate-like Growth 

We have applied the algorithms derived above to examine a key problem in snow 
crystal growth (and, by extension, other examples of faceted crystal growth) - 
understanding how small changes in extrinsic parameters like temperature and 
supersaturation can produce rather dramatic changes in the resulting crystal 
morphologies. Using our new LCA growth algorithms, we discovered an inter- 
esting example of a morphological transition in the growth of plate-like snow 
crystals. 

We modeled the growth of small snow crystal plates in the limit of cylindrical 
symmetry, so we could use the 2D algorithm in cylindrical coordinates described 
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above. This is much simpler computationally than the full 3D problem, and the 
approximation is a reasonable one for small plates without dendritic branching 
j!8j . Instead of six prism facets, the cylindrical model produces one continuous 
"facet" which is the perimeter of the plate. Thus, for example, hollowing of the 
six prism faces is replaced by hollowing of the perimeter " facet" . Aside from 
geometrical factors of order unity, we believe the cylindrical model is a good 
approximation for the growth of simple plate-like snow crystals. 

In our model we used a grid of N r — 200 by N z = 100 pixels with A£ = 1, 
corresponding to a physical grid spacing of Ax = 0.15 /mi. The boundary 
condition on the z — plane guaranteed symmetry about that plane. We began 
each growth run with a single ice pixel at r = z = 0, and we ran the code until 
the maximum crystal radius reached 20 fim. For the basal surface, we assumed 
that the growth was nucleation-limited with a(a) = min[l,Aexp(— cto/ct)] with 
A = 2 and ct — 0.021, following the latest ice crystal growth measurements 
[16] . Accurate prism growth measurements do not yet exist, so we assumed 
nucleation-limited growth with A = 5 and cto = 0.01. For non-facet surfaces 
we chose an = 0.7 for the kink site with B = 3 and a = 1 for all sites with 
B > 3. We believe that the morphological transition we observed is insensitive to 
the exact values of these growth parameters, as long as they produce plate-like 
structures with nucleation-limited growth on both facets. 

Our results are shown in Figures 2 and 3 as we varied only the supersatura- 
tion CToo in the model. The figures show a clear morphological transition from 
simple, thin-plate growth at low CToo to concave (hollowed) growth at interme- 
diate CToo and convex growth at high CToo. 

The growth of thin plates at low CToo (number 1 in the figures) is easy to 
understand from simple considerations of the growth dynamics. At low (Too, a 
is small and the growth is largely limited by attachment kinetics, so ct sa Coo at 
the crystal surface. The relative growth velocities of the prism and basal faces 
is then 

Vprism siprism r/ \ / l { ac\\ 

~ —. exp [(ao,basal — &0,prism) /CToo J (4U) 

Vbasal -^-basal 

Since ao^asai > o'o.prism for our case, this ratio increases rapidly as CToo de- 
creases, producing thinner plates at lower CToo- 

As CToo increases, the higher ct values at the corners causes hollowing of both 
the prism and basal facets, as seen in Figure 3, number 2. For this concave 
growth, we see that steps are generated at the edges of the basal facets and 
propagate inward as the crystal grows. 

As CToo increases still further, the growth undergoes a rather sudden transi- 
tion to convex growth, seen in Figure 3, number 3. This transition is accom- 
panied by a reduction in crystal thickness and growth time, as seen in Figure 
2. These two are related in that a thinner crystal requires less mass for a given 
radius, and thus grows to a given radius (20 /im in this case) in less time. For 
this convex growth, steps are generated at the centers of the basal facets and 
propagate outward, in contrast to the concave case. 

At the highest CToo shown, the basal growth has increased until additional 
structure is seen at the center of the basal facet. Similar structures appear when 
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the growth rates of the prism and basal facets are comparable, and this growth 
behavior reflects the fact that ocbasai / otprism increases with increasing a^ as 
the surface values of a begin to exceed the critical supersaturations cfq for the 
facets. 

The transition from concave to convex growth as the supersaturation in- 
creases is well-known in snow crystal growth pQ , because it is one of the most 
studied examples of faceted crystal growth. Other crystal systems may ex- 
hibit similar morphological transitions as the growth drive is increased. The 
morphologies that result certainly depend on the attachment kinetics, so the 
detailed modeling must be tailored to each individual system. 

The LCA algorithms described above have proven to be quite robust and nu- 
merically stable in these calculations, while being simple to code with reasonably 
fast execution times. It is straightforward to change the attachment kinetics and 
other growth parameters to investigate different regimes. The LCA method is 
particular well-suited to modeling the horizontal propagation of macrosteps on 
a flat surface, which appears to be a key element of faceted growth. Further 
studies are needed to identify whatever numerical idiosyncrasies are inherent in 
the LCA method (such as the intrinsic anisotropy described above), but so far 
it appears to be a powerful and flexible approach for modeling faceted crystal 
growth. 
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Figure 1: A mapping of a hexagonal grid onto a rectangular grid. Corresponding 
pixels are numbered, showing the arrangement of nearest neighbors in each case. 
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Figure 2: Modeled growth time (top panel) and final thickness (lower panel) of 
a crystal as a function of the supersaturation at the outer boundary. For all 
crystals the growth was terminated when the diameter reached 40 /an. Growth 
parameters are described in the text. The numbers in the lower panel correspond 
to those in Figure 3 showing the crystal morphology. 
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Figure 3: Full cross-sections of four plate-like crystals at the end of their growth, 
at which time the diameter was 40 fira. The numbers beside each crystal corre- 
spond to those in Figure 2. This shows the transition from simple plate growth 
(1) to concave growth (2) to convex growth (3 and 4) as (Too is increased. 
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